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ABSTRACT 

We examine several different simplified approaches for modelling the chemistry of CO 
in three-dimensional numerical simulations of turbulent molecular clouds. We compare 
the different models both by looking at the behaviour of integrated quantities such 
as the mean CO fraction or the cloud-averaged CO-to-H2 conversion factor, and also 
by studying the detailed distribution of CO as a function of gas density and visual 
extinction. In addition, we examine the extent to which the density and temperature 
distributions depend on our choice of chemical model. 

We find that all of the models predict the same density PDF and also agree very 
well on the form of the temperature PDF for temperatures T > 30 K, although at 
lower temperatures, some difference become apparent. All of the models also predict 
the same CO-to-H2 conversion factor, to within a factor of a few. However, when we 
look more closely at the details of the CO distribution, we find larger differences. 
The more complex models tend to produce less CO and more atomic carbon than the 
simpler models, suggesting that the C/CO ratio may be a useful observational tool 
for determining which model best fits the observational data. Nevertheless, the fact 
that these chemical differences do not appear to have a strong effect on the density 
or temperature distributions of the gas suggests that the dynamical behaviour of the 
molecular clouds on large scales is not particularly sensitive to how accurately the 
small-scale chemistry is modelled. 

Key vifords: galaxies: ISM - ISM: clouds - ISM: molecules - molecular processes 



1 INTRODUCTION 

Carbon monoxide (CO) is a key constituent of the gas mak- 
ing up the giant molecular clouds (GMCs) that are the site 
of almost all Galactic star formation. Although almost all 
of the molecular gas mass within a CMC is in the form of 
molecular hydrogen (H2), this material is extremely difficult 
to observe directly, as the temperatures within a CMC are 
too low to excite even the lowest accessible rotational transi- 
tion of the H2 molecule, the quadrupole transition between 
the J = 2 and J = rotational levels. On the other hand, 
CO does become rotationally excited at CMC temperatures, 
owing to the much smaller energy separation between its ro- 
tational levels. It is readily observed in the millimeter, and is 
arguably the single most important observational tracer of 
the state of the gas within the clouds. It also plays an impor- 
tant role as the main molecular coolant of the gas over a wide 



range in densities (see e.g. Neufeld, Lepp fc Melnick||1995 1. 
It is therefore important to understand the distribution of 
CO within GMCs, and how this relates to the underlying 
gas distribution. 



Unfortunately, this is not a simple task. It has been un- 
derstood for a long time that the gas within GMCs is typ- 
ically not in chemical equilibrium (see e.g. Leung, Herbst 



& Huebner 19841. More recently, numerical modelling has 



highlighted the fact that the CO abundance within a given 
parcel of gas within a GMC is a complex function of the 
density and temperature of the gas, the local radiation field 
(modulated by absorption by gas elsewhere in the cloud), 



and possibly also the dynamical history of the gas (Glover 



et al.l2010). It has also become increasingly clear that molec- 



ular clouds are dynamically complex objects, dominated by 



supersonic, turbulent motions (see e.g. Mac Low & Klessen 
|2004[ for a comprehensive overview) , and that we ignore the 
effects of this dynamical complexity at our peril. Therefore, 
if we want to be able to accurately model the formation and 
destruction of CO within a GMC, the use of detailed models 
that treat both the chemistry and the turbulent dynamics 
in an appropriately coupled fashion seems unavoidable. 

However, this presents us with a serious numerical chal- 
lenge. Even if we focus only on the chemistry of hydrogen, 
helium, carbon and oxygen, and ignore other elements such 
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as nitrogen or sulphur, there are still a very large number 
of possible reactions and reactants that we could include 
in our chemical model. For example, in the 2006 release of 



the UMIST Database for Astrochemistry (Woodall et al 



2007), there are 2115 such reactions, involving a total of 



172 different chemical species. In particular, it is the large 
number of chemical species that must potentially be han- 
dled that is the main source of our problems. If we model 
the chemical evolution of the gas by using a set of ordi- 
nary differential equations (ODEs) to represent the rates of 
change of the chemical abundances, then we find that the 
required set of ODEs is generally extremely stiff, containing 
processes with a wide range of different intrinsic timescales. 
To evolve this set of ODEs in a numerically stable fashion, 
we must therefore either evolve them explicitly in time using 
extremely small timesteps, which is impractical, or evolve 
them implicitly in time. However, if we adopt an implicit 
approach, then we find in general that the computational 
cost of solving the implicit ODEs scales as the cube of the 
number of species involved, A'^pec- Therefore, when Nupcc is 
large, the cost of solving the chemical rate equations can eas- 
ily come to dominate the total computational cost, putting 
strict limits on the size of the problems that can be tackled. 
In practice, values of Ngp^c as large as 14 can be handled 
in three-dimensional turbulence simulations at the present 
day (see e.g. Glover et al.|2010 note that while their model 
includes 32 distinct chemical species, all but 14 of these are 
assumed to be in chemical equilibrium, or are trivially deriv- 
able from conservation laws), but even in this case, solving 
the chemical rate equations can take up as much as 90% of 
the total runtime of the simulation. Scaling up from here 
to the Napec = 172 subset of the UMIST database men- 
tioned above remains computationally impractical in three- 
dimensional simulations. 

For this reason, a number of authors have sought to 
reduce Napec to a more tractable value by identifying and 
retaining only the most important of the reactions involved 
in the formation and destruction of CO, as well as mak- 
ing other simplifications that are discussed in later sec- 
tions. When identifying the key reactions, there is an obvi- 
ous trade-off between complexity (and hence computational 
cost) and completeness: the larger the number of possible 
chemical pathways that we include, the slower the code will 
run. There are models available in the literature that span a 
range of different complexities, but until now there has been 
no comparison of the results produced by these different ap- 
proaches. 

In this paper, we present results from a detailed compar- 
ison of several different approximate methods for modelling 
CO formation and destruction in turbulent molecular clouds 
that span a range of different complexities. We are interested 
in particular in establishing which results from simulations 
of turbulent clouds are highly sensitive to the choice of sim- 
plified chemical model, and which are relatively insensitive. 



2 SIMULATION SETUP 

The basic setup of the simulations used here is similar to 



that described in Glover et al. ( 2010 1 and Glover & Mac Low 
( |2011[ ). We use the ZEUS-MP code to model driven magne- 
tohydrodynamical turbulence in a cubical volume, with side 



length L = 20pc, using periodic boundary conditions for the 
gas. The turbulence is driven using the algorithm described 
in Mac Low et al. (19981, and the amplitude of the driving 



is chosen in order to maintain the rms turbulent velocity 
at Virns = 5 km s~^. The initial magnetic field strength is 
-Bo = 5.85 /iC, and the field is assumed to be initially uni- 
form and oriented parallel to the 2;-axis of the simulation. 
The effects of self-gravity are not included. The simulations 
are run for t = 1.8 x lO^'^s ~ 5.7Myr, or approximately three 
turbulent eddy turnover times. We showed in |Glover et al.| 
(20101 that this is long enough to allow the CO distribution 



to settle into a statistical steady state. 

We perform two sets of simulations (hereafter denoted 
as set 1 and set 2) with different initial densities: no = 
100 cm""^ for set 1 and no — 300 cm~^ for set 2, where 
no is the initial number density of hydrogen nuclei. The cor- 
responding mean extinctions are therefore Av = 3.3 and 
Ay = 9.9, respectively, where we have assumed a conversion 
factor from hydrogen column density to dust visual extinc- 
tion of 5.348 X 10~^^ mag cm^, appropriate to dust in the 
cold interstellar medium (ISM) . We therefore expect that in 
both cases, detectable amounts of CO will be produced in 
the clouds, even though the mean mass- weighted CO abun- 
dance should differ by more than an order of magnitude 
between the two sets of runs ( |Glover fc Mac Low|20lT| ). 

In both sets of simulations, we adopt elemental abun- 
dances (by number, relative to hydrogen) of 2;hc = 0.1, 
xc ~ 1.41 X 10~^ and xq = 3.16 x 10~* for helium, carbon 
and oxygen, respectively. We assume that the hydrogen, he- 
lium and oxygen are initially atomic, and that the carbon 
is initially in singly ionized form throughout the simulation 
volume. 

We set the initial gas temperature to 60 K and the ini- 
tial dust temperature to 10 K in both sets of simulations. To 
model the subsequent heating and cooling of the gas and the 
dust, we use a modified version of the [Glover et aL\ ([2010} 
thermal model, as described in that paper and in the Ap- 
pendix. Set 1 and set 2 consist of six simulations each, one 
for each of the chemical models described in Section [S] In 
the models that include the effects of cosmic rays, we adopt 
a rate ("h = 10~^^s~'^ for the cosmic ray ionization of atomic 
hydrogen. Finally, to model the effects of the external ultra- 
violet radiation field, we assume that its strength and shape 



are the same as those of the standard Draine ( 1978 1 field. We 



illuminate each side of the simulation box using the unat- 
tenuated field, and model attenuation within the box due 
to dust absorption, H2 self-shielding, CO self-shielding and 
the shielding of CO by H2 using the 'six-ray' approximation 



introduced by Nelson & Langer (19971. Full details of this 



procedure can be found in Glover et al. (20101 



3 CHEMICAL MODELS 

3.1 Glover et al. (2010) [GlOg, GlOng] 

The first model we consider is a slightly modified version 
of the [Glover et al.| ( [2010[ ) chemical model. This consists 
of 218 reactions amongst 32 species, and so although it is 
the most complicated of the approximate models included 
in this study, it nevertheless still represents a considerable 
simplification compared to a model including the full set of 
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reactions from e.g. the UMIST Database for Astrochemistry. 



dn, 



CO 



However, Glover et al. (2010 1 demonstrated that this simpli- 



fied model could accurately reproduce the results obtained 
with a far more comprehensive model for several ID test 
problems, while at the same time being small enough to be 
usable in a three-dimensional simulation. 

We have made one significant change to the |Glover et al.| 



(20101 chemical network: the inclusion of an optional treat- 
ment of the effects of the recombination of H^, He^, C"*", 
and ions on the surfaces of charged dust grains. We 



treat these processes using the formalism of [Weingartner fc] 
Draine (2001 1, which includes the effects of very small dust 



grains and poly cyclic aromatic hydrocarbons. These grain 
surface processes are not included in any of the other mod- 
els that we examine, and in order to understand what effect 
they have on the outcome of the simulations, we have per- 
formed runs both with and without them. In the remain- 
der of the paper, we denote these simulations as GlOg and 
GlOng, respectively. 

We have also made a number of modifications to our 
treatment of the thermodynamic behaviour of the gas, in 
order to improve our ability to model very cold gas. These 
improvements are used for all of the models examined here 
and full details of them are given in the Appendix. 

3.2 Nelson &: Langer (1997) [NL97] 

A much simpler approximation is given by the model pro- 



posed by Nelson & Langer ( 19971, which they used to study 
the dynamics of low-mass (M = 100-400 Mq) molecular 



clouds. In their study. Nelson & Langer ( 1997 1 assume that 



all of the hydrogen is already in the form of H2, and focus 
their attention on the conversion of singly ionized carbon, 
C^, to carbon monoxide, CO. Their approximation involves 
the assumption of direct conversion from C"^ to CO, and vice 
versa, which allows them to ignore any intermediate species 
(such as neutral atomic carbon, C). They assume that the 
conversion of C^ to CO is initiated by the formation of an 
intermediate hydrocarbon radical (e.g. CH or CH2), which 
they denote as CH^,. This may then react with oxygen to 
form CO, or be photodissociated by the interstellar radia- 
tion field. Once CO has formed, it is then only destroyed by 
photodissociation, yielding C and O, but the neutral carbon 
produced in this way is assumed to be instantly photoion- 
ized, yielding C"*". Since the formation of the hydrocarbon 
radical will typically involve a slow radiative association re- 
action as the initial step, such as the formation of CH^ via 

C++R2^ CH+ + 7, (1) 

[Nelson fc Langer] ( |1997[ ) assume that this is the rate limiting 
step for the formation of CO, and write the rate equation 
for the CO number density asQ 



^ In actual fact, Nelson & Langer give slightly different forms for 
Equations|2]and[4] writing n in place of njjj in Equation[2] (equa- 
tion 18 in their paper) and vice versa in Equation |4] (equation 20 
in their paper). However, this appears to be a typographical er- 
ror, as one can see by considering the behaviour of the system 
when the H2 number density and the UV field strength are both 
zero. According to the logic of the NL97 model, the formation 
rate of CO in this case should also be zero, since H2 is required in 
order to form the intermediate hydrocarbon radical. However, if 



dt 



— koUQ+n^^P — Tconco, 



(2) 



where nc+ is the number density of C^ ions and is the 
number density of hydrogen molecules. In Equation [2] fco 
is the rate coefficient for the formation of the intermediate 
CHa; ion or radical, which Nelson fc Langer ( 1997[ l give as 
fco = 5 X 10~^^ cm^ s~^, and Fco is the photodissociation 
rate of CO, given in their model as 

Leo = 10~'°G'oexp(-2.5^v) s~\ 



(3) 



where Go is the strength of the ultraviolet radiation field 



in units of the Habing ( 1968 1 field. In our simulations, we 



adopt the |Draine| ( 1978 1 parameterization of the interstellar 
ultraviolet radiation field and so Go = 1.7. The variable /3 
in Equation [2] represents the proportion of the GRx that 
successfully forms CO. This is given in the [Nelson fc Langer| 
(19971) model by 



kixo 



kiXQ + FcH^/n' 



(4) 



where n is the number density of hydrogen nuclei, xq is 
the fractional abundance of atomic oxygen, fci is the rate 
coefficient for the formation of CO from O + CHa;, given by 
Nelson & Langer as fci = 5 x 10~^° cm'' s~^, and FcHx is the 
photodissociation rate of CHa, , given by 



TcH, = 5 X 10""'Goexp(-2.5ylv) s" 



(5) 



We have implemented this treatment of the carbon 
chemistry, with a couple of minor changes. In place of the 
rate assumed by [Nelson fc Langer[ ( [1997[ ) for Fco, we use 
a rate Fco = 2 x 10~^^Go exp(— 2.5Av)/sh s"^, where fsh 
is a shielding factor that quantifies the effects of CO self- 
shielding and the shielding of CO by H2 Lyman-Werner 
band absorption. This is the same rate coefficient as that 
used in the GlOg model, and is based on work by [van[ 
Dishoeck fc Black[ (IT9881) and [Lee et al.[ ((19961). We have 



made this substitution in an effort to minimize any differ- 
ences between the models that arise purely from differences 
in the rate coefficients adopted. The influence of rate coef- 
ficient uncertainties on molecular cloud chemistry has been 
studied in detail elsewhere (see e.g. [ Millar et al.| 1988 Wake^ 
lam, Herbst fc Sclsis 2006 ; Wakclam ct al. 2010) and is not 



our primary focus here, as we arc interested more in the in- 
fluence of the design of the chemical network itself. Since we 
have adopted a larger value for Fco, we have also adopted 
a larger value for Fch^ , so as to keep the ratio of FcHx/Fco 
in the optically thin limit the same as in [Nelson fc Langer[ 
(T997|. This yields a value for FqHx that is more in keeping 
with the rates adopted for CH and CH2 photodissociation 



and photoionization in the GlOg model (see Glover et al 



2010 Appendix A). To compute the shielding factor f^h, we 



use our standard six-ray treatment to determine the H2 and 
CO column densities, and then convert these to a shielding 



factor using the data tabulated in Lee et al. ( 1996 1 



In addition, and unlike Nelson fc Langer 



(1997), we do 



not assume that the hydrogen is completely molecular, but 



one uses 



the original expressions given in [NeIson fc Langer[ | [1997l l 
for and /3, one finds that they do not show this behaviour 

— instead, /3 — > 1 in the limit that njjj and the predicted CO 
formation rate remains larger than zero. 
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instead follow the evolution of the H2 and H"'" abundances 
explicitly using the same hydrogen chemistry network as 
Glover & Mac Low ( 2007a|b \ . The abundance of neutral 



atomic hydrogen, H, then follows from a simple conservation 
law. We include the effects of H2 self-shielding and dust 
shielding using the same six-ray treatment as in the GlOg 
model. 



3.3 Nelson & Langer (1999) [NL99] 

In a later paper, Nelson & Langer suggested an alternative 



approximation for modelling the formation of CO (Nelson 



& Langer 19991. This was designed for a similar purpose 



as the Nelson & Langer ( 1997 1 approximation, but is con- 



siderably more sophisticated. Notably, it allows for the for- 
mation of CO via multiple pathways. In addition to the for- 
mation channel involving the composite hydrocarbon radical 
CHj; (which should be understood to represent both CH and 
CH2) that forms the basis of the NL97 model, the Nelson 
fc Langer] ( |1999[ ) model (hereafter NL99) also allows for CO 
formation via the composite oxygen species OH^, (represent- 
ing the species OH, H2O, O2 and their ions) as well as via 
the recombination of HCO"'' . In addition, photodissociation 
is no longer the only fate for the CO: the network also in- 
cludes the conversion of CO to HCO"^ by proton transfer 
from Hj^, and its destruction by dissociative charge transfer 
from ionized helium: 



CO + He+ ^ C+ + O -f He. 



(6) 



A further notable difference between the NL97 and NL99 
models is the fact that the latter model tracks the abun- 
dance of neutral atomic carbon, rather than just C"*" and 



CO. Finally, Nelson & Langer ( 1999 1 also include in their 



model a small number of reactions involving a species they 
denote as M that represents the combined effects of low ion- 
ization potential metals such as Mg, Fe, Ca and Na, which 
become the dominant atomic charge carriers in very shielded 
regions of the cloud. The full list of reactions included in the 
Nelson & Langer ( 1999 1 model is given in Table [1] 

Many of the reactions in the NL99 model are also in- 
cluded in the GlOg model, and for these reactions we adopt 
the same rate coefficients in both models, for the reasons dis- 
cussed previously. For the reactions not in the GlOg model - 
notably, those involving CHa;, OHa; or M, we adopt the same 
rate coefficients as in NL99. For the elemental abundance of 
M, which is not included in any of our other chemical mod- 



els, we adopt the same value as in Nelson & Langer ( 1999 1, 
i.e. XM.tot = 10~^, and we assume that it is initially fully 
ionized. 

We also supplement the list of reactions given in Table[l] 
with those used in the |Glover &; Mac Low| f2007a|b| ) network 
for hydrogen chemistry, and once again include the effects 
of shielding using our standard six-ray approach. As in the 
case of the NL97 model, we include the effects of CO self- 
shielding and the shielding of CO by H2 in order to allow us 
to make a fair comparison with the GlOg model. 



3.4 Keto & Casein (2008) [KCOSe, KCOSn] 

The final two approximations that we consider in this study 



Table 1. Reactions in the NL99 cliemical model 
Reaction Notes 



H2 



H+ H2 ^ H+ H 
He -I- c.r. He+ 



C- 
O- 



H 



CO- 
He+ 
He+ 
C+- 
C+- 
O-hCHx 
C-hOHx 
He+ -I- e- 
H+ +e- 
C+ +e- 
HCO+ + 
M+ +e- 
- M - 



•CHx 
' OHx 

h;^ HCO+ 



H 



H 

e~ 
H2 
H2 

H2 

H2 ^ He -I- H -I- H+ 
- CO ^ C+ -I- O -I- He 
H2 ^ CHx H 
OHx HCO+ 
CO + H 
CO + H 
' He -I- 7 
H2 +H 
C + 7 
^ CO- 
M + 7 
H+ + M -> M+ + H2 - 
C -I- 7 -> C+ -(- e- 

CHx + ^^c + n 

CO -I- 7 c -f o 
OHx + 7 ^ O + H 
M-I-7 -> M+ -|-e~ 
HCO+ + 7 -s> CO + H+ 



Notes: 1: These two reactions are combined into a single pseudo- 
reaction in NL99, as it is assumed that all of the H^ formed by 
the first reaction is immediately consumed by the second. 2: M 
represents the combined contributions of several low ionization 
potential metals, such as Mg, Fe, Ca and Na. 

veloped for the study of the thermal balance in dense prestel- 



H 



lar cores. As in Nelson & Langer (19971, they assume that 



CO forms primarily via an intermediate hydrocarbon radi- 
cal, explicitly assumed in this case to be CH2, and that the 
formation of this radical is the rate limiting step in the for- 
mation of CO. Unlike [Nelson fc Langer] ( |1997| ), they do not 
account for photodissociation of the CH2, and so write the 
rate equation for the CO number density as 



dnr 



dt 



fcRAWH2'T'C+ ~ LcOWCO, 



(7) 



where /cra is the rate coefficient describing the formation 
of CH2 from C''' and H2 (via radiative association to form 
CHJ, which is then rapidly converted to CH2) and Fco is 
the photodissociation rate of CO, the only destruction mech- 
anism for CO that is included in their models. Unlike iNel-l 



son & Langer (19971, they do not assume that the carbon 



produced by CO photodissociation will be instantly pho- 
toionized, and hence write the rate equation for the neutral 
carbon number density as 



dnc 



rco?ico - Tcnc, 



(8) 



where Fc is the photoionization rate of atomic carbon. The 
rate equation for the C^ number density then follows as 

dnc+ 



dt 



— Fcnc — fcRA?iH2?ic+- 



(9) 



In their study, Keto & Caselli ( 2008 1 assume chemical equr 



are based on the work of Keto & Caselli ( 2008 1, and were de- 



librium, and hence eliminate all of the time derivatives from 
the above set of equations, yielding a set of coupled algebraic 
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equations. If one also makes use of the conservation equation 
relating the fractional abundance of C"'', C and CO to the 
total fractional abundance of carbon, a;c,tot, 



Table 2. Computational performance of the various approaches 



X(^+ + XC + xco = a;c,tot, 



(10) 



then these algebraic equations are simple to solve for the 
equilibrium fractional abundances of C^, C and CO. One 
obtains: 



Xc+ 

a;c,tot 

a:c,tot 
and 

xco 



1+^ + 



Xc+ 



Xc+ 



Xc 
Xc+ 



Xc+ Xc+ J 



a;c,tot \Xq+ j y Xq+ Xq+ ^ 
where the ratio of CO to C"'' is given bjj^ 

xco _ fcRAflHa 

Xc+ Tco 

and the ratio of C to C^ is given by 



xc 

Xc+ 



(11) 



(12) 



(13) 



(14) 



(15) 



Keto & Caselli (20081 adopt numerical values for fcaA, Tco 



and Fc from Tielens & Hollenbach (19851, but as in the 



other models we examine, we instead use values from |Glover| 
et al. (20101 in order to minimize any differences arising 



purely from differences in the adopted rate coefficients. As 
in the other models, we account for the shielding of CO by 
CO and H2 when computing Fco- 

We have implemented this equilibrium treatment of the 
carbon chemistry and coupled it with our standard treat- 
ment of the non-equilibrium hydrogen chemistry, described 
in Glover & Mac Low (2007a I. We denote this implementa- 



tion as KC08e (for "equilibrium") in the remainder of the 
paper. We have also implemented a non-equilibrium carbon 
chemistry based on the same set of reactions, but which uses 
the time-dependant rate equations (Eqs. [7||9]l as its basis, 
rather than the equilibrium abundances. This implementa- 
tion is denoted below as KC08n (for "non-equilibrium"); 
note also that Keto & Caselli (20101 use a similar non- 
equilibrium scheme. 



4 PERFORMANCE 

We begin our comparison of these various approaches to 
modelling CO formation and destruction in turbulent gas 
by examining their relative computational performance. In 
Table [2] we list the computational time required (in units 



2 Note that Equation 5 in [Keto fc Caselli| l |2008[ l, which gives 
xco/^C+ rco/Fc is incorrect. The authors have kindly con- 
firmed to us that this is due to a typographical error: T'co/^C 
actually corresponds to the ratio of the C and CO abundances, 
xc/^COy rather than xco/^c+ printed. The authors have ver- 
ified that this error only appears in the journal article, and not 
in the original numerical modeling, and hence it does not afliect 



Method 


Approximate runtime (CPU hours) 




Run 1 


Run 2 


GlOg 


1190 


1360 


GlOng 


880 


1240 


NL97 


110 


140 


NL99 


310 


460 


KCOSe 


100 


130 


KCOSn 


120 


150 



of CPU hours) to perform the simulations discussed in this 
paper. All of the simulations were run on 32 cores on kolob, 
an Intel Xeon Quad- Core cluster at the University of Hei- 
delberg]^ Note that as we made no special efforts to ensure 
that the computational workload of the cluster was com- 
pletely identical during each run, these numbers should be 
treated with a certain amount of caution - they are rea- 
sonably indicative of the computational performance of the 
different approaches, but could easily be uncertain at the 10- 
20% level. Nevertheless, some clear trends are apparent. The 
three approaches that attempt to model the carbon chem- 
istry with only one or two additional rate equations (NL97, 
KCOSe and KC08n) are clearly the fastest, but do not differ 
significantly amongst themselves in terms of required run- 
time. This likely just reflects the fact that the time taken 
up by the chemistry in these simulations does not dominate 
the total computational cost. 

The NL99 model is approximately a factor of three 
slower than these simpler models, reflecting its significantly 
greater complexity, but is a factor of three faster than the 
even more complex GlOg and GlOng models. The slowdown 
as we go from NL99 to GlOg or GlOng is roughly in line 
with what we expect, given an A'^'^pcc scaling for the cost 
of the chemistry: A^spoc = 10 in the NL99 model (nine 
non-equilibrium species, plus the internal energy), while 
A^'apcc = 15 in GlOg and GlOng (14 non-equilibrium species, 
plus the internal energy), and so the expected slowdown is 
a factor of 15^/10^ ~ 3.4. On the other hand, the slowdown 
between e.g. NL97 and NL99 is much smaller than we might 
expect, given that Nspec = 4 in NL97, consistent with the 
chemistry not being the dominant cost in the simplest mod- 
els. Finally, the difference in runtimes between models GlOg 
and GlOng is likely due to the fact that the grain-surface 
recombination rate coefficients are somewhat costly to cal- 



culate, at least if one uses the Weingartner & Draine (2001 \ 



any of the results presented in Keto & Caselli I 2008 I 



fitting functions, and their calculation has not yet been sig- 
nificantly optimized within the current version of the code. 

In terms of the limits that the requirement of following 
the chemistry places on the size of simulation that can be 
performed, it is worth bearing in mind that a factor of two 
increase in spatial resolution in a three-dimensional ideal 
MHD Eulerian simulation will generally lead to a factor of 
sixteen increase in runtime: the number of resolution ele- 
ments increases by a factor of eight, while the Courant con- 
dition causes the maximum timestep to decrease by a factor 
of two, so that twice as many timesteps are required to reach 



For full technical details, see |http://kolob.ziti.uni" 
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the same physical time. In simulations where the chemistry 
is the dominant computational cost, and is also subcycled 
(i.e. evolved on a timestep smaller than the magnetohydro- 
dynamical timestep), then this last factor of two increase 
can often be avoided, since the change in spatial resolution 
does not directly affect the number of chemical substeps that 
must be taken. However, even in this case, the bottom line 
is that a factor of two improvement in spatial resolution will 
lead to an order of magnitude increase in runtime. There- 
fore, the difference in performance between the simplest and 
most complex chemical networks is roughly equivalent to a 
factor of two in spatial resolution: a simulation performed 
with 256"^ zones using model GlOg will take roughly as long 
as a simulation performed with 512^ zones using one of mod- 
els NL97, KC08e or KC08n. 



5 RESULTS 

5.1 Chemical abundances: time evolution 

Having discussed the relative computational performance of 
the various approaches, and shown that, as expected, the 
simpler models are considerably faster than the more com- 
plex ones, we now look in more detail at the behaviour they 
predict for the chemical abundances. We begin with some of 
the simplest quantities that we might expect the models to 
be able to reproduce, the total masses of atomic carbon and 
CO formed in the simulation as a function of time. A con- 
venient way to express this is in terms of the mass-weighted 
mean abundances of these species. We can define the mass- 
weighted mean abundance of a species p as 



J2i i k ^p(*' ^)^(*' J' k)AV{i,j, k) 



(16) 



where Xp{i,j,k) is the fractional abundance (by number, 
relative to the number of hydrogen nuclei) of species m, 
p{i,j,k) is the mass density in zone {i,j,k), AV{i,j,k) is 
the volume of zone {i,j,k), Mtot is the total mass of gas 
present in the simulation, and where we sum over all grid 
zones. It is simple to convert from {a;p)M to Mp (the to- 
tal mass of species p in the simulation) using the following 
equation: 



Mp = 



-Mtot{a;p)M, 



(17) 



where mn is the mass of a hydrogen atom and muc is the 
mass of a helium atom. 

In Figure [l^, we plot the evolution with time of the 
mass- weighted mean CO abundance, {a;co)M, in the six sim- 
ulations in set 1, our low density simulations. It is immedi- 
ately obvious from this figure that the different methods do 
not agree on the amount of CO that is formed in the gas. 
Although there is good agreement between closely related 
models (e.g. KC08n and KCOSe, or GlOg and GlOng), there 
is roughly an order of magnitude difference in the CO frac- 
tions predicted by the full set of models. Furthermore, it 
is clear that the three simplest models (NL97, KCOSe and 
KC08n) find a qualitatively different evolution for the mean 
CO abundance that the three more complex models (NL99, 
GlOg and GlOng), predicting a more rapid rise in {a;co)M 
and a significantly larger final value. 




I I 
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Figure 1. (a) Time evolution of the mass- weighted mean CO 
abundance in the simulations in set 1, representing the evolution 
of a low density cloud with mean hydrogen number density no = 
lOOcm"^. (b) As (a), but for the simulations in set 2, which model 
the evolution of a higher density cloud, with no = 300 cm~^. 




Figure 2. (a) Time evolution of the mass-weighted mean abun- 
dance of atomic carbon in the simulations in set 1, the low density 
cloud. Note that no line is plotted for the NL97 model, as this 
model does not track the abundance of atomic carbon, (b) As (a), 
but for the simulations in set 2. 
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The behaviour of {a;co)M in our higher density runs in 
set 2, illustrated in Figure [TJj, further supports our finding 
that we can divide our models into two subsets that predict 
qualitatively different behaviour for the growth of the CO 
mass fraction. Runs NL97, KC08e and KC08n form almost 
all of their CO at t < 1 Myr, and thereafter predict that 
(a;co)M should barely evolve. On the other hand, models 
NL99, GlOg and GlOng find that CO forms over a more 
extended period, and predict that {a::co)M reaches a statis- 
tical steady state only at t > 5 Myr. Runs NL97, KC08e 
and KC08n agree very well on the final amount of CO pro- 
duced, predicting values for (a:co)M that differ by only a 
few percent. Runs NL99, GlOg and GlOng also agree well 
at times t > 1 Myr on the amount of CO produced, but 
predict a final value for {a;co)M that is only 60% of the size 
of that predicted by runs NL97, KC08e and KC08n. The 
only difference between runs GlOg and GlOng is the inclu- 
sion of grain-surface recombination in the former, and the 
good agreement between the CO distributions produced by 
these two models therefore indicates that the inclusion of 
this process has only a minor effect on the CO abundance. 
The good agreement between the results of these two models 
and the NL99 model is more of a surprise, and suggests that 
the composite species CHj, and OH^; in the NL99 model do 
a good job of reproducing the behaviour of the intermediate 
hydrocarbons and oxygen-bearing species that are tracked 
directly in models GlOg and GlOng. 

In Figure [2j we examine the evolution of the mass- 
weighted mean abundance of atomic carbon, {a;c)M, in our 
simulations. In this case, we compare only five models, as 
the NL97 chemical model does not contain atomic carbon, 
and so cannot make any prediction whatsoever about its 
abundance. It is immediately apparent from the Figure that 
there is a far larger degree of disagreement for the C abun- 
dances than for the CO abundances, particularly for our low 
density simulations. In these, we see three different evolu- 
tionary tracks. Run GlOg forms a large amount of atomic 
carbon very quickly, reaching a value of {a;c)M ~ 2 x 10~^ 
within only a few thousand years, an interval so short that 
it is not well represented in the Figure. Following this, the 
atomic carbon abundance continues to increase for the rest 
of the run, but at a much, much slower pace, reaching a final 
abundance (a;c)M ~ 5.5 x 10~^ by the end of the simula- 
tion. In runs GlOng and NL99, there is also a rapid growth in 
the atomic carbon abundance at very early times, but in this 
case, this phase of rapid growth stops once {a;c)M ~ 2x 10~®. 
From this point until t ~ 1 Myr, the atomic carbon abun- 
dance continues to increase significantly, albeit at a much 
slower pace than at the very beginning of the simulation, 
while for t > IMyr, there is little further evolution of {a::c)M. 
Finally, in runs KC08e and KC08n, the atomic carbon abun- 
dance does not display the very rapid growth at t <^ 1 Myr 
seen in the other runs. Although there is significant growth 
in the atomic carbon abundance at t < 1 Myr, the character- 
istic timescale is a significant fraction of a megayear, rather 
than the few thousand years found in the other models, and 
the amount of atomic carbon formed is much smaller. At 
t > 1 Myr, there is a clear change in behaviour; the growth 
in {xc)m slows significantly, and it appears to have reached 
a steady-state by the end of the simulation. The final abun- 
dance of atomic carbon is considerably smaller than in the 
other runs: it is roughly a factor of thirty smaller than in 



than in runs GlOng or NL99, and more than a factor of 100 
smaller than in run GlOg. 

In our higher density simulations, the size of the dis- 
agreement between run GlOg and runs GlOng and NL99 is 
somewhat smaller, although significant differences remain. 
As in the lower density case, we see a very rapid growth in 
the atomic carbon abundance at very early times in these 
models, followed in this case by a slow decline as the atomic 
carbon is incorporated into CO molecules. The KC08e model 
also predicts a rapid rise in {2;c)m followed by a slow decline, 
but the amount of atomic carbon produced in this model is a 
factor of 100 or so smaller than in the other models. Finally, 
model KC08n predicts a somewhat slower rise in (a;c)M at 
early times, with a characteristic timescale of about 1 Myr, 
followed by a decline in the atomic carbon fraction at times 
t > 1 Myr that is almost indistinguishable from that in the 
KC08e model. 

These discrepancies in the behaviour of the atomic car- 
bon abundance in the various models are actually rela- 
tively easy to understand. In the KC08n and KC08e mod- 
els, atomic carbon is produced only as an outcome of CO 
photodissociation, rather than directly from C"*" by recom- 
bination. Since the formation rate of CO is relatively slow, 
relying as it does on a radiative association reaction, this 
means that the equilibrium abundance of atomic carbon is 
very small when the visual extinction is small. Therefore, in 
the KC08e model, the growth of the carbon abundance is 
regulated by the appearance of regions with high visual ex- 
tinctions. In run 1, regions with a large enough visual extinc- 
tion to allow for the production of non-negligible amounts 
of C and CO are created by the turbulent restructuring of 
the gas on a timescale of roughly 1 Myr. In run 2, on the 
other hand, the higher mean density allows some regions to 
have a significant visual extinction (and hence a high equi- 
librium abundance of atomic carbon) even at t — 0. In the 
KC08n model, the same consideration applies with regards 
to the dependence on visual extinction, but in addition, the 
growth of the atomic carbon abundance is also regulated by 
the time taken to form CO, which itself is of the order of 
1 Myr or longer in much of the gas. 

In the other three models, the most effective way to 
convert C*" to C is by direct recombination: gas-phase re- 
combination in models GlOng and NL99, and a mix of 
gas-phase and grain-surface recombination in model GlOg. 
The recombination time of fully ionized carbon in gas with 
n — 100 cm~'^ and T = 60 K is roughly 0.2 Myr, assuming 
that the ionized carbon itself is the dominant source of the 
required electrons, and since the equilibrium C/C^ ratio at 
t = is much smaller than unity, it requires only a small 
fraction of a recombination time for the C/C^ ratio to reach 
its equilibrium value. Subsequent changes in the atomic car- 
bon abundance are driven by two main effects: the restruc- 
turing of the gas by the turbulence, which puts more of the 
gas mass into dense, high Ay regions where the equilibrium 
C/C^ ratio is larger, and the conversion of C into CO in this 
same dense, well-shielded gas. This latter process is respon- 
sible for the fall-off in the atomic carbon abundance in the 
high density versions of runs GlOg, GlOng and NL99 that 
occurs at t > 1 Myr (see Figure [2]). The larger C/C^ ratio 
that we find in run GlOg compared to runs GlOng and NL99 
is an obvious consequence of the inclusion of grain-surface 
recombination, which is particularly effective when the ratio 
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of the UV field strength Go to the mean density no is small, 
as it is in our simulations. 

Finally, we note that although it would be simple to 
investigate the evolution of the total mass of C"'" in a similar 
fashion to our analysis of C and CO above, we have chosen 
to omit this, as in practice, C*", C and CO between them 
contain almost all of the available carbon in models GlOg, 
GlOng and NL99 (and must contain 100% of it in the other 
models), and so there is nothing new to be learned from 
studying the time evolution of the C^. 



5.2 Chemical abundances: dependence on density 
and visual extinction 

We can gain more insight into the differences between the 
various models tested here by looking in more detail at their 
predictions for the distribution of CO at the end of the sim- 
ulations. We know from previous work ( [Glover et aT||2010[ ) 
that regardless of whether we examine them as a function 
of density or visual extinction, we always find considerable 
scatter in the CO abundances in a turbulent cloud. The rea- 
son for this is that the CO abundance is sensitive to both 
density and visual extinction, but these two quantities are 
only poorly correlated within turbulent clouds. Nevertheless, 
despite this scatter it can still be useful to examine averaged 
quantities, such as the total fraction of carbon represented 
by C^, C or CO, and how this varies as a function of density 
or visual extinction. 

Therefore, in Figure [Sj we show how the fraction of 
carbon in the form of C"'' , C or CO varies with density in the 
simulations. Let us focus initially on Figure |3^, which shows 
the results from the simulations in set 1. At densities below 
the mean density of 100 cm~'^, all of the models predict that 
almost all of the carbon will be C"*", and hence agree very well 
in their predictions of the C^ abundance. Ifowever, when we 
look at the abundances of C and CO at these low densities, 
we find substantial disagreement between the results of the 
different models. Let us start by considering the abundance 
of atomic carbon. This is entirely absent in the NL97 model. 



but represents 0.1% of the total carbon at n = 100 cm in 
the KC08 models, roughly 1% in models GlOng and NL99, 
and closer to 10% in model GlOg. The differences in the 
predictions for the atomic carbon abundance coming from 
various models persist over a wide range of densities, and 
at high densities, where atomic carbon is abundant, they 
also cause significant differences in the C^ fraction, which, 
for instance, falls off more rapidly with increasing density in 
run GlOg than in the other runs. Turning to the CO fraction, 
we find that runs GlOg, GlOng and NL99 agree well with 
each other on the behaviour of the CO fraction, as do runs 
NL97, KC08e and KC08n, but that there is a significant 
difference in the behaviour of these two sets of runs. All 
of the models predict a sharp rise in the CO fraction with 
increasing number density, and agree that the fraction of 
carbon in CO should be close to 100% for number densities 
of order lO^cm""^. However, at densities below lO'^cm"'', we 
find that the CO fraction in runs NL97, KCOSe and KCOSn 
is systematically larger than in the other three runs. For 
example, at n = 1000 cm~^, the three simple models predict 
a CO fraction of roughly 80-90%, while runs GlOg, GlOng, 
and NL99 predict a value that is closer to 10%. 

In Figure [Sb, which shows the behaviour of the simu- 



lations in set 2, we find very similar results. In this case, 
the transition from C^ to C to CO takes place over a wider 
range of densities, and the gas becomes CO dominated at a 
lower density than before. This is a consequence of the higher 
mean extinction of the gas - in these simulations, there is 
more gas at all densities that is well shielded from the UV 
background and hence can maintain a high CO abundance. 
Nevertheless, the qualitative behaviour of the CO fraction 
as a function of density remains the same as in our lower 
density runs. Models NL97, KC08e and KCOSn continue to 
agree well over the whole range of densities plotted, as do 
models GlOg, GlOng and NL99, with the first set of models 
predicting systematically higher CO fractions than the sec- 
ond set. At n > 5000 cm"'', all six models agree that the gas 
should be completely CO-dominated. 

In Figure |4] we show how the fraction of carbon in the 
form of C"*", C or CO varies as a function of the effective vi- 
sual extinction, defined for any given cell in our simulations 
as 



A 



V,cff = 



2.5 



In 




g-2-5^v(^P+) ^ g-2-5Av{a:p_) 



(18) 



where Av{xp+) is the visual extinction of material between 
that cell and the edge of the volume in the positive direction 
along the and so forth. The choice of the factor of 

2.5 occurs because in our models, the CO photodissociation 
rate scales with the visual extinction Ay as exp(— 2.5j4v). 
The value of j4v,off defined in this fashion corresponds to the 
visual extinction used in our code, in the context of our six- 
ray approximation, for computing the CO photodissociation 
rate. 

Figure |4] displays a number of familiar features. Models 
GlOg, GlOng and NL99 all agree well on the evolution of 
the CO fraction with increasing visual extinction, just as 
they did on its evolution with increasing density. Models 
GlOng and NL99 also agree on the behaviour of the C and 
C^ fractions at all but the highest visual extinctions, while 
model GlOg predicts a much larger atomic carbon fraction 
at ylv.off ~ 1-2 than in the other two models, and hence a 
correspondingly smaller C^ fraction. Models NL97, KC08e 
and KCOSn agree well with each other regarding the CO 
fraction, but produce more CO at low visual extinctions than 
the other three runs. 

Figure |4j3 also displays a curious feature in the plot 
of the atomic carbon fraction at very low Ay^^s in runs 
GlOg, GlOng and NL99. Below an effective visual extinc- 
tion of roughly 0.4, the atomic carbon fraction increases with 
decreasing Ay.eff in these three runs. However, further in- 
vestigation shows that this feature is a numerical artifact 
related to our use of the six-ray approximation. As previ- 
ously discussed in [Glover et al.| ( [2010[ ) , the effective visual 
extinction of material right at the edge of the computational 
volume is higher than it should be, owing to the poor angu- 
lar sampling of the radiation field, and in practice only few 
zones very close to the edge have effective visual extinctions 
Ay^aS < 0.5. Because of this, when we compute the mass 
fractions for Ay^^s < 0.5, we are averaging over only a small 
number of zones, and hence are far more sensitive to the 
effects of outliers than we would be at higher Ay.eff. In this 
particular case, the anomalous behaviour of the atomic car- 
bon fraction is due to the contribution from a small clump of 




Figure 3. (a) Fraction of the total available carbon found in the 
form of C+ (dashed lines), C (dash-dotted lines), or CO (solid 
lines) in the simulations in set 1 , plotted as a function of density. 
Note that as model NL97 does not include atomic carbon, no line 
is plotted in this case, (b) As (a), but for the simulations in set 
2. 



dense gas located at the edge of the simulation volume. The 
higher density of this gas enables recombination to be 
more effective, and allows it to have a higher atomic carbon 
abundance than the rest of the gas at this visual extinction. 
The effect is more pronounced in run GlOg because of the 
inclusion of grain-surface recombination in that model. 



5.3 CO column densities and integrated 
intensities 

Although informative, the quantities that we have examined 
so far have the significant disadvantage that they are not 
observable in real molecular clouds. It is therefore useful to 
examine whether the differences between the models that we 
have already discussed will lead to clear differences in the 
observable properties of the CO distribution. 

We begin by examining the behaviour of the CO col- 
umn density distribution. CO column densities are directly 
observable in the ISM along low extinction {Av ~ 2 or less) 
sightlines to bright background stars, where they can be ac- 
curately determined by UV absorption measurements (see 
e.g. [Sheffer et al.|[2007[ ). For higher extinctions, this tech- 
nique is no longer effective, and one typically relies on CO 
emission, which is an accurate tracer of the CO column den- 
sity only for transitions which are optically thin and uni- 



formly excited (se e e.g. Pineda, Caselli & Goodman 2008 



Shetty et al. 2011 1. Since much of the '^CO emission coming 



from Galactic GMCs is optically thick, information on the 
CO column density distributions in these clouds generally 
comes from observations of rarer CO isotopologues such as 
"CO or C^®0, which are optically thin over a much wider 
range of cloud column densities. 



Figure 4. (a) As Figure|3^, but showing the variation in the 0+, 
C and CO fractions as a function of visual extinction ylv for the 
simulations in set 1. (b) As (a), but for the simulations in set 2. 



The probability density functions (PDFs) for the CO 
column density in our two sets of runs are plotted in Fig- 
ure [5] In our low density runs, we see that the PDFs are 
relatively flat, indicative of there being a roughly equal prob- 
ability of selecting any particular CO column density within 
a relatively wide range. This is quite unlike the PDF of to- 
tal column density, which has the characteristic log-normal 
shape ubiquitously found in simulations of supersonic turbu- 
lence (see e.g. Ostriker, Stone fc Gammie|2001 1, and demon- 
strates that the CO is not a particularly accurate tracer of 
the underlying density distribution in our low density simu- 
lations (see also Shetty et al.|20lT for more on this point). 
In runs NL97, KC08e and KCOSn, we begin to see the influ- 
ence of this underlying structure between CO column den- 
sities of lO^'' cm~^ and a few times lO'^^ cm~^, where there 
is a clear peak in the PDF, but no such feature is visible in 
the CO column density distributions produced in the other 
three runs. In our higher density runs, on the other hand, 
the influence of the underlying density structure of the gas 
is far more pronounced, with all of the runs showing a clear 
peak in the PDF at high CO column densities. The PDFs 
can best be understood as the superposition of two different 
features: a log-normal portion at high A^'co, corresponding 
to lines of sight along which most of the carbon is in molecu- 
lar form, leading to a CO column density that simply traces 
the total column density, and an extended tail at much lower 
A'^co that corresponds to lines of sight along which the CO 
is significantly photodissociated. 

Comparing the predictions of the different models, we 
see that once again the runs can be separated into two 
distinct sets. Runs GlOg, GlOng and NL99 agree reason- 
ably well with each other, particularly for the higher den- 
sity cloud, but produce results that differ significantly from 
those found in runs NL97, KCOSe and KCOSn. The latter 
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Figure 5. (a) CO column density PDFs for the runs in set 1. (b) 
As (a), but for the runs in set 2. 



models produce narrower CO column density PDFs, with a 
more pronounced peak at high Ai'co, and with much smaller 
probabilities at low A'^co- One consequence of this is that the 
CO is a better tracer of the underlying density distribution 
in runs NL97, KC08e and KC08n than in runs GlOg, GlOng 
and NL99. 

We have also computed estimates of the frequency- 
integrated intensities of the J = 1 — > line of CO cor- 
responding to these CO column densities, using the same 



technique as in Glover & Mac Low (2011 1. To briefly sum- 



marize, we assume that the CO is in local thermodynamic 
equilibrium (LTE), and is isothermal, with a temperature 
equal to a weighted mean temperature for the gas, computed 
using the CO number density as the weighting function: 



11, CO ~ 



(19) 



where we sum over all grid zones i,j,k. We also assume that 
the CO linewidth is uniform, and is given by Av = 3kms~^ 
(which is roughly equal to the one-dimensional velocity dis- 
persion we would expect to find in a gas with an RMS tur- 
bulent velocity of 5 kms~^). Given these assumptions, we 
can relate the CO column density to the optical depth in 
the J = 1 line using ( |Tielens||2005| ) 



Aioc^gi 
8tvuIq go 



fo 



exp 



iVco 
Av ' 



(20) 



where Aio is the spontaneous radiative transition rate for the 
J = 1 — >■ transition, vio is the frequency of the transition, 
iJio — hvio is the corresponding energy, go and gi are the 
statistical weights of the J = and J = 1 levels, respectively, 
and fo is the fractional level population of the J = level. 
We can then convert from no to the integrated intensity, 
Wco, using the same curve of growth analysis as in [Pineda^ 
Caselh & Goodmanl (|2008D: 



Figure 6. (a) PDF of the integrated intensity in the J = 1 — > 
line of 12 CO in the runs in set 1. The integrated intensities are 
estimates, calculated using the technique described in |GIover &| 
[Mac Low] l |201l"] | and summarized in Section |5.3[ (b) As (a), but 
for the runs in set 2. 



Wco = TbAi 



2/3(f)df, 



(21) 



where Tb is the brightness temperature of the line, which 
we assume is simply equal to T'mcan,co, and P is the photon 
escape probability, given by 



[l-exp(-2.34r)]/4.68r r sC 7, 
4r [in (r/y^)] 



1/2 



r > 7. 



(22) 



Although the estimates of Wco generated by this procedure 
are probably accurate to within only a factor of a few, this 
level of accuracy is sufficient for our purpose here. 

In Figure [6] we plot PDFs of integrated intensity for 
the simulations. At Wco < 30Kkms~^, these PDFs closely 
resemble the CO column density PDFs, while at higher inte- 
grated intensities, we see a clear peak in the PDF even when 
there is no corresponding feature in the CO column density 
distribution at high A*'co- This behaviour can be easily un- 
derstood as being due to the change from having optically 
thin CO emission at Wco ~ 30K km s"^ and below, to hav- 
ing optically thick emission at higher Wco - In the optically 
thin regime, Wco oc A^co, and so we see similar behaviour 
in both PDFs. In the optically thick regime, on the other 
hand, Wco increases only slowly with increasing A'^co, and 
so we find a 'pile-up' of values at these intensities (see also 
the more detailed discussion of this effect in |Shetty et al.| 
20TT1). 



5.4 The C/CO ratio 

Another potentially useful quantity for distinguishing be- 
tween different models is the ratio of the column density of 
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Table 3. Comparison of C/CO column density ratios 



Table 4. Comparison of the CO-to-H2 conversion factor, Xqq, 
at the end of the runs 



Method 



Nc/Nco 





Run 1 


Run 2 


GlOg 


9.3 


0.68 


GlOng 


2.2 


0.28 


NL99 


1.5 


0.23 


KCOSe 


0.01 


0.002 


KC08n 


0.01 


0.002 



atomic carbon to that of CO. A number of groups have at- 
tempted to measure this value (see e.g. Ingalls et al.|[l997 



Ikeda et al.||2002 Bensch et al.|2003 l, typically finding val 



ues in the range of 0.1 to 3, with the higher values coming 
from lower column density translucent clouds, and the lower 
values coming from higher column density dark clouds. Al- 
though a comprehensive comparison between our results and 
these observational determinations is outside of the scope of 
our current paper, a simple comparison nevertheless proves 
illuminating. 

In Table [S] we list the ratio of the mean column density 
of atomic carbon, Nc, to the mean column density of CO, 
Nco, that we obtain for each of our runs, with the exception 
of the two NL97 models, which do not track atomic carbon. 
We see that two of the models - GlOng and NL99 - produce 
values for the C/CO column density ratio that are broadly 
in line with the values observed in real clouds. On the other 
hand, model GlOg in run 1 produces a significantly higher 
value than is observed, suggesting that our treatment of C""" 
recombination on grains may actually overestimate the rate 
at which this process occurs in the real ISM. In this con- 



text, it is interesting to note a recent study by Liszt (2011 1 
that also finds indications that we do not currently under- 
stand the role that grain-surface recombination of C^ plays 
in the transition from C^ to C to CO. Finally, we see from 
Table [3] that the KCOSe and KCOSn models produce signif- 
icantly smaller C/CO ratios than are seen in real clouds, 
likely because these models neglect the effect of gas-phase 
recombination of C"*". 



5.5 The CO-to-H2 conversion factor 

Before we conclude our study of the details of the CO dis- 
tribution in our simulations, it is interesting to examine 
whether the differences in the amount of CO formed in the 
various runs, and in its spatial distribution, lead to signif- 
icant differences in the CO-to-H2 conversion factor, Xco- 
This is conventionally defined as the ratio of the H2 col- 
umn density to the integrated intensity of the J — 1 ^ 
transition of CO, i.e. 



Xco 



Wco' 



(23) 



Observations of Galactic GMCs show that Xco appears to 
be roughly constant from cloud to cloud, with a mean value 
of Xco = 2 X 10^°cm~^ K'^ km~^ s (see e.g. [Dame et aT 



2001), and in Glover & Mac Low (20111 we showed that 



our turbulent cloud models reproduce this behaviour, pro- 
vided that the mean visual extinction of the clouds exceeds 
a threshold value of Av ~ 2-3. 



Method 


Xco [10^° cm- 


2 (Kkms-i)-i] 




Run 1 


Run 2 


GlOg 


2.67 


1.30 


GlOng 


2.01 


1.23 


NL97 


0.57 


0.90 


NL99 


1.47 


1.17 


KC08e 


0.50 


0.89 


KC08n 


0.49 


0.87 



To estimate Xco for the simulations considered in this 
paper, we start by finding the mean of the distribution of 
integrated intensities that we have already computed for 
each simulation. Armed with a mean integrated intensity 
for each simulation, we next compute the mean H2 column 
density for each simulation, following which we can obtain 
our estimate of Xco simply by taking the ratio of these two 
mean quantities. The values obtained using this procedure 
are listed in Table |4l 

In the low density case, models NL97, KCOSe and 
KCOSn yield values for Xco that are a factor of a few smaller 
than the standard Galactic value, while the other three runs 
yield values that are in much better agreement with the ob- 
servations. However, given that our estimates of Xco here 
are likely only accurate to within a factor of a few, it is 
unclear whether the difference between predictions of the 
NL97, KCOSe and KCOSn models and the observations is 
really meaningful. In any case, in the high density case we 
find much better agreement between the models, with all of 
them now producing values that are somewhat smaller than 
the Galactic value. The reason for this similarity is that all 
of the runs agree reasonably well (to within a factor of two) 
on the peak value of the integrated CO intensity found in the 
simulation, with the significant differences in the Wco dis- 
tribution occurring for values considerably below the peak. 
Our estimate of Xco is dominated by the contribution of 
lines of sight with integrated intensities close to the peak 
value, and is insensitive to the behaviour of lines of sight 
with low Wco, and so this diagnostic tells us little about 
the low intensity sightlines. 



5.6 Gas temperature 

C+, C and CO are all important coolants at the tempera- 
tures and densities found within molecular clouds, and so 
differences in the predicted distributions of these species 
may have an effect on the thermal evolution of the gas and 
on its temperature structure at the end of the simulation. 

We have investigated this by examining the time evo- 
lution of the mass-weighted mean gas temperature, plotted 
in Figure [7] In the low density runs, the mean temperature 
drops very rapidly from our initial value of 60 K to a value 
of approximately 30 K. Thereafter, it evolves very little over 
the course of the simulation. The six different chemical mod- 
els predict very similar values for the mean temperature. The 
most discrepant outcome is for model GlOg, which predicts 
a final mean temperature that is roughly 1 K cooler than 
in the other runs. This result is easy to understand given 



12 



Glover & Clark 



, , 1 1 1 1 

-a) — — 


^ ' 1 ' ' '- - 




GIO ~ 




GlUng " 




NL97 ; 




NL99 ~ 




KC08e - 




1 KCOBn - 


~: — 1 1 1 1 — 


H 1 1 1 \ — : 



b) 




_L 



2 4 
Time [Myr] 




Temperature [K] 



Figure 7. (a) Time evolution of the mass-weighted mean gas 
temperature in the simulations in set 1. (b) As (a), but for the 
simulations in set 2. 



Figure 8. (a) Mass-weighted temperature PDF for run 1. (b) As 
(a), but for run 2. Note that in both cases, the results for the 
NL97 and KCOSe runs are barely distinguishable from those for 
the KC08n run. 



the results for the mean mass- weighted abundances of C , 
C and CO in these runs, discussed earher in Section |5.1| 

is the dominant form of carbon in all of these runs, and 
hence will be the dominant coolant. Differences in the abun- 
dances of C or CO will therefore have only a minor impact 
on the energy balance of the gas, and hence will have little 
effect on the mean gas temperature. However, in the case of 
model GlOg, the atomic carbon abundance becomes large 
enough that its contribution to the cooling of the gas can 
no longer be neglected. Since the energy separation of the 
fine structure levels in neutral carbon is significantly smaller 
than in singly ionized carbon, the former is a more effective 
low temperature coolant than the latter, and so increasing 
the C abundance at the expense of the C^ abundance leads 
to an overall lower temperature for the gas. 

In the higher density runs, the mean temperature again 
drops very rapidly initially, before stabilizing at a value of 
around 17-19 K. In this case, the difference between the 
models is more pronounced, with runs NL97, KC08 and 
KC08n predicting the highest temperatures, and run GlOg 
predicting the lowest. However, at the end of the simulations, 
the difference between the runs is only 2 K, or roughly 10%. 
The differences can again be easily understood on the basis 
of our prior results for the mean abundances. In this case, 
the C and CO abundances are much higher, and so C^ is 
no longer the dominant radiative coolant. Instead, the cool- 
ing is typically dominated by CO (runs NL97, KCOSe and 
KCOSn) or by a mixture of atomic carbon and CO (the other 
three runs). There is a clear inverse correlation between the 
amount of atomic carbon present in the gas and the final 
mean temperature. 

We have also examined the temperature distribution of 
the gas at the end of the simulations, as shown in Figure [Sj 
where we plot mass-weighted PDFs of the gas temperature 



for the various runs. In both the low and the high density 
cases, the PDF at T > 30 K shows almost no dependence 
on the chemical model, consistent with it representing gas 
which is dominated by C^ in all of the models. Below 30 K, 
however, clear differences become apparent. In the low den- 
sity case, runs GlOng, NL99, KCOSn, KCOSe and NL97 all 
show a clear peak at a temperature of around 20 K, but 
in run GlOg, this peak occurs at the lower temperature of 
10 K. Runs NL97, KCOSe and KCOSn all produce very little 
gas cooler than 10 K, but in the other three runs we find 
a significant fraction of gas with temperatures in the range 
of 5-10 K. In the high density run, the behaviour changes 
slightly. Runs NL97, KCOSe and KCOSn still have a clear 
feature in the PDF at T ~ 20 K, but the true peak is now 
found at a temperature slightly larger than 10 K. In runs 
GlOg, GlOng and NL99, on the other hand, the peak is at 
the slightly lower temperature of 9 K, with this feature being 
more pronounced in run GlOg than in the other two runs. 
The behaviour of these three runs at temperatures below 
this peak value is very similar. 

The differences between the low temperature behaviour 
in the various runs can be relatively easily understood. In 
the low density runs, little CO is formed, and so in most 
of the gas, the dominant cooling mechanisms are C^ and C 
fine structure emission, with some contribution from dust 
in the densest regions. The energy separation of the 
ground state and the ^P3/2 excited state of C"*" is approxi- 
mately 92 K, and hence it cannot easily cool the gas below 
about 15-20 K, owing to the exponential suppression of the 
cooling rate. On the other hand, atomic carbon has a separa- 
tion of only 24 K between its two lowest lying energy levels, 
and so remains an effective coolant down to temperatures 
of order 5 K. Therefore, chemical models producing greater 
quantities of atomic carbon will tend to produce lower gas 
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temperatures. For instance, as we have already seen, run 
GlOg produces more atomic carbon in intermediate density 
gas (n ~ 1000 cm"^) than runs GlOng or NL99. In the low 
density case, little CO is formed, as we have already seen, 
and so most of the 10 K gas produced in this case comes from 
regions that are dominated by cooling from neutral carbon. 
Consequently, run GlOg produces significantly more of this 
10 K gas than runs GlOng or NL99. On the other hand, run 
NL97 produces no neutral carbon, and hence only very little 
10 K gas. In the higher density simulations, which have much 
higher mean CO abundances, much more of the 10 K gas is 
situated in regions dominated by CO cooling, and so the 
differences between runs GlOg, GlOng and NL99 are much 
smaller than in the low density case. The lack of neutral 
carbon in run NL97, and the fact that it underproduces CO 
at low densities mean that it still disagrees with these three 
models, and produces less very cold gas. The behaviour of 
runs KCOSe and KC08n is very similar to that of run NL97 
because the small amount of neutral carbon produced in 
runs KCOSe and KCOSn is not enough to significantly affect 
the thermal balance of the gas, and as we have already seen, 
the CO and C^ abundances in these runs are very similar 
to those in run NL97. 



5.7 Gas density 

Finally, we have investigated whether the differences in the 
temperature distributions examined above lead to any sig- 
nificant differences in the density distribution of the gas, 
as quantified by the density PDF. This is of particular im- 
portance for understanding the level of accuracy in the mod- 
elling of the CO chemistry that is necessary in order to allow 
us to model star formation accurately, as there are a number 
of theoretical models that argue that it is the shape of the 
density PDF that determines the stellar initial mass func- 
tion (jPadoanT^ordlm 

20021 iHennebelle & Chabrier»2008! "20091) or the star forma- 



tion rate ( [Krumholz fc Mc Kco 2005| [Padoan fc Nordlun d 
|2011[ ) in a molecular cloud. 

In Figure |9^, we plot the density PDFs for the simula- 
tions in set 1. In Figure |9]3, we show a similar plot for the 
simulations in set 2. We see that in both cases, over most of 
the density range spanned by the PDF, there is essentially 
no difference between the results of any of the runs, indicat- 
ing that at these densities, the physical structure of the gas 
is insensitive to the CO abundance. Minor differences start 
to become apparent in the high density tail of the PDFs, 
but even here the differences are small. In addition, they 
occur at densities that are not well resolved in our current 
simulations, and we cannot be completely certain that they 
would persist in higher resolution simulations. 



Padoan, Nordlund & Jones ( 1997 1 presented results in- 
dicative of a relationship between the dispersion as of the 
logarithmic density contrast s = ln(p/po) (where po is the 
mean gas density), and the volume- weighted RMS Mach 
number M, finding that 

o-f = In (l b'M^) , (24) 
with b ~ 0.5. Recent work by |Federrath, Klessen fc Schmidt] 



oT 




10 100 10^ 10" 

n [cm"^] 



1 2008 1 has shown that the proportionality parameter b is 
sensitive to the relative strength of the solenoidal and com- 
pressive modes in the forcing field used to drive the turbu- 



Figure 9. (a) Mass-weighted density PDF for the runs in set 
run 1. Minor differences are apparent in the high density tail, at 
n > 5000 cm~^, but over the rest of the density range, all six runs 
agree very well, (b) As (a), but for the runs in set 2. Again, there 
is very good agreement between the runs at almost all densities, 
with only a few minor differences becoming apparent in the high 
density tail. 



lence, and ranges from 6 ~ 1/3 for purely solenoidal forcing 
to fe ~ 1 for purely compressive forcing. In each case, the 
effects of the gas temperature enter through the dependence 
of CTs on the volume-weighted RMS Mach number. In our 
simulations, the dominant contributions to A4 come from 
low density, warm gas that is dominated by C"'' cooling. 
The differences in the temperature distribution of the cold, 
denser material have little influence on M. We find that A4 
varies by no more than 1-2% from run to run in both the 
low density and high density cases, thereby explaining why 
we see so little variation in the density PDF from run to 
run. 



6 DISCUSSION 

6.1 Understanding the differences in CO 
production 

Our comparisons in the previous section have demonstrated 
that the GlOg, GlOng and NL99 models agree very well re- 
garding the time evolution of the CO abundance and the 
spatial distribution of the CO. The most significant differ- 
ence between the treatment of the carbon chemistry in these 
three models turns out to be the inclusion of grain-surface 
recombination in model GlOg, or more specifically, the re- 
combination of C^ ions on grain surfaces. The inclusion of 
this process has a significant effect on the ratio of C"*" to C 
in the gas, and therefore indirectly affects the temperature 
distribution, since C is a much more effective coolant than 
C^ at gas temperatures below 20 K. 
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The KCOSe, KC08n and NL97 models also agree well 
with each other, but disagree with the GlOg, GlOng and 
NL99 models in two major respects. First, they predict a 
shorter formation timescale for the CO, particularly in the 
high density run. Second, they predict systematically larger 
values for the CO fraction as a function of density or vi- 
sual extinction. Both of these disagreements between the 
two sets of models result from the same underlying cause. In 
the KCOSe, KCOSn and NL97 models, the rate-limiting step 
in the formation of CO is the initial radiative association 
between C""" and H2, a reaction that has a rate coefficient 
fco = 5 X 10~^^ cm'^ s~^ in the NL97 model, and a similar 
value in the other two models. The rate per unit volume 
at which CHj ions (or CJix radicals) form in this model is 
therefore given by 



-Ri = 5 X 10 riQ+TiHa- 



(25) 



In the GlOg, GlOng and NL99 models, on the other hand, 
there are multiple routes leading to the formation of CO and 
hence the situation is somewhat more complicated. In prac- 
tice, the main contributions come from two main reactions 
pathways, one initiated by the radiative association of C"*" 
with H2, as above, and a second initiated by the reaction of 
atomic carbon with Hg". The latter reaction occurs rapidly, 
and the rate-limiting step for this second reaction pathway 
is hence not the reaction between C and Hj^, but rather the 
formation of the ion itself. This occurs as a consequence 
of the cosmic ray ionization of H2, which yields ions that 
then rapidly react with H2 to form HJ: 



H+ -I-H2 



H++H. 



(26) 



In regions where at least a few percent of the carbon is 
ionized, the free electron abundance is roughly equal to the 
abundance of ionized carbon, and the dominant destruction 
mechanism for the ions is dissociative recombination: 



H++e 



H2+H, 



(27) 
(28) 



In these conditions, the net rate at which CHj ions are 
produced by reactions between C and is therefore given 
by 



R2 = 2CHnH2 



CH+ nc 



(29) 



fcdr nQ+ ' 

where is the cosmic ray ionization rate of atomic hydro- 
gen, fcdr is the rate coefficient for the dissociative recombi 



nation of H, ions and k, 



cm 



reaction 



c + m 



CH.T + H. 



is the rate coefficient for the 



(30) 



In cold gas, k^-^^+ /fcdr — 0.01, and hence 

7?2 ^ 2 X 10-''Ch.i7 nn^ — , (31) 

"-0+ 

where Ch,i7 = Ch/10^^^s~^. (Recall that in our simulations, 
^H,i7 = 1). The net rate of formation of CHJ ions in models 
GlOg, GlOng and NL99 is approximately equal to Ri -\- R2, 
since other processes make only minor contributions, while 
in models NL97, KCOSe and KCOSn, the rate of formation 
of CHJ ions is simply Ri. Note, however, that since _Ri 
depends on 71^+ , which is generally larger in models NL97, 



KCOSe and KCOSn than in the other models, the size of 
the rate _Ri in the simple models is generally larger than in 
the more complex models. To avoid confusion, it is useful 
to distinguish between these two cases by writing the rate 
in the simple models (NL97, KCOSe, KCOSn) as Ri,s and 
writing the rate in the more complex models (NL99, GlOng, 
GlOg) as Ri,c. 

We next consider the conditions in which the rate of 
CHJ formation in the complex models, given by Ri^c + R2, 
is bigger than the rate in the simple models, given by iii,s. 
Starting with the inequality, 

Ri,c + R2 > Ri,s, (32) 
we can subtract iii^c from both sides, giving us 



R2 > (-Ri.s — ^l,c 



(33) 



If we denote the number density of C'*' in the simple mod- 
els as nc+ g, and denote the same quantity in the complex 
models as nQ+ ^, then we can write this inequality as 



2xlO-'^CH,i7nH2 



nc 



> 5 X 10 ^®nH,(n, 



■c+ 



c)(34) 



where we have assumed that the H2 abundance is the same 
in both cases. If the CO abundance is small compared to 
the abundances of C^ and/or C, then in the simple models, 
the number density of C^ ions will be roughly equal to the 
number density of carbon atoms in all forms, i.e. nQ+ ^ ~ 
wctot- In the complex models, on the other hand, we have 
n^+.c + nc — nc,tot, and hence in these conditions 71^+ ^ — 
7ic+,c — ^c- We can use this fact to simplify Equation 34 
which reduces to the following constraint on nQ+ 



T,c+,c<4xl0 Ch,17. 



(35) 



In practice, this inequality is rarely satisfied in the gas in 
our simulations. For example, in solar metallicity gas with 
n = 300cm~'^, it is satisfied only when the ratio of atomic to 
ionized carbon is of the order of a thousand to one. There- 
fore, in general, the rate at which CHj ions are produced in 
the more complex models is slower than the rate at which 
they are produced in the simpler models. 

This fact explains most of the difference that we see be- 
tween our two sets of chemical models. When nQ+ ^ nc, 
all of the models produce CHj ions (and hence also CO 
molecules) at a very similar rate. In models GlOg, GlOng 
and NL99, however, the physical conditions that allow large 
CO fractions to be produced (high density and high visual 
extinction) also allow the recombination of C''' to be far 
more effective than the photoionization of atomic carbon, 
meaning that the equilibrium ratio of C to C^ in these re- 
gions strongly favours atomic carbon. As the carbon recom- 
bines (which occurs on a timescale that is short compared 
to the CO formation timescale), the CO formation rate de- 
creases for the reasons outlined above. The same effect does 
not occur in the NL97, KCOSe or KCOSn models because 
none of these models include the effects of C^ recombina- 
tion: model NL97 ignores atomic carbon entirely, while in 
models KCOSe and KCOSn, it is produced only by the pho- 
todissociation of CO. Therefore, the CO formation rate in 
models NL97, KCOSe or KCOSn is in general larger than in 
models GlOg, GlOng and NL99, explaining the shorter CO 
formation timescales and larger CO fractions we find in the 
former models with respect to the latter. 
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r CONCLUSIONS 

In this paper, we have examined several different approaches 
to modeUing the chemistry of CO in three-dimensional nu- 
merical simulations of turbulent molecular clouds. The sim- 
plest approaches (NL97, KCOSn, KC08e) include only a 
very small number of reactions and reactants, with the 
KC08e model making the additional simplifying assumption 
of chemical equilibrium. The NL99 adds a number of ad- 
ditional reactions and reactants, but still remains relatively 
simple, owing to its use of artifical species (CHa,, OHa,) as 
stand-ins for a much wider range of real molecular ions, radi- 
cals and molecules. Finally, the most complex models exam- 
ined here (GlOg and GlOng) add a few more non-equilibrium 
reactants, a large number of additional reactants that are as- 
sumed to be in chemical equilibrium, and a large number of 
additional reactions. The complexity of these two models 
lies close to the upper limit of what is practical to include 
in a three-dimensional simulation at present, but still repre- 
sents a significant simplification when compared with state- 
of-the-art one-zone or one-dimensional chemical models (see 
e.g. the PDR codes discussed in |Rollig et al.|2007 l. 

As one would expect, the simpler one makes the chem- 
ical model, the faster it becomes to compute the chemical 
evolution of the gas. Models NL97, KC08e and KCOSn al- 
low the simulations to run roughly an order of magnitude 
faster than the reference model, while model NL99 yields a 
runtime three times longer than the simple models, but still 
three to four times shorter than the models based on lGloverl 
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Looking in more detail at the results of the simple net- 
works, we find that they tend to produce more CO than 
the more complex networks, for the reasons outlined in the 
previous section. It is also clear that one area in which the 
simpler models fall down is in their treatment of atomic 
carbon. This is ignored entirely in the NL97 model, and is 
significantly underproduced in the KCOSe and KCOSn mod- 
els, owing to the omission of the effects of C"^ recombination 
from these models. It should however be noted that the lKetol 
|fc Caselli] ( |2008[ ) chemical model was designed for use in a 
study of the thermal balance of prestellar cores with mean 
densities that are more than an order of magnitude larger 
than the mean densities considered in our simulations. In 
these conditions, the gas is CO-dominated, atomic carbon 
is not an important coolant, and the'Keto & Caselli (2008) 
model does a good job of capturing the essential physics of 
the gas. It performs more poorly in our study only because 
we are applying it to a range of physical conditions outside 
of those for which it was designed. 

Turning our attention to the more complex networks, 
we find that models NL99 and GlOng produce very similar 
results in all of our comparisons, and that the differences 
between the results of these models and that of model GlOg 
are easily understandable as consequences of the inclusion 
of grain-surface recombination of C"*" in the latter, which 
significantly affects the relative abundances of C and C''' 
produced in the gas. We also note that the inclusion of this 
process leads to a significant overproduction of atomic car- 
bon and hence yields values for the C/CO ratio that are 
significantly higher than those observed in real molecular 
clouds. Models NL99 and GlOng, on the other hand, yield 



values for the C/CO ratio that appear to be consistent with 
observations. 

Given that model NL99 produces very similar values to 
model GlOng for the C and CO abundances, but in only 
one-third of the computational time, we suggest that the 
former is a better choice than the latter if one is primarily 
interested in the C^ to C to CO transition, or in the de- 
tails of the resulting CO distribution. On the other hand, 
the more complex GlOng model is a better choice if one is 
also interested in the distributions of species not tracked in 
the NL99 model, such as CH, CH+, OH, H2O or O2, and is 
willing to pay the additional computational cost for tracking 
these species. Furthermore, if all one is interested in is an 
approximate picture of which regions are likely to have high 
CO abundances (as may be of interest in a large-scale simu- 
lation of the Galactic ISM, such as Dobbs et al.|2008 l, then 
any of these models are likely to do a reasonable job, as the 
differences between the CO distributions that they predict 
are not enormous and are likely smaller than the uncertain- 
ties introduced by the inevitably limited resolution possible 
within large-scale simulations. For this kind of modelling, 
it would therefore make sense to choose one of the simpler, 
faster models, such as NL97, KCOSe or KCOSn. 

Finally, we have examined the effects of the choice of 
chemical network on the density and temperature distribu- 
tions in our model clouds. We find that the influence of the 
chemistry is surprisingly limited. In both the low and the 
high density model clouds, much of the cloud volume is filled 
with warm gas (T > 30 K). Most of the carbon in this warm 
material is in the form of C*", and hence C"*" cooling domi- 
nates. All of the chemical networks that we examined do a 
good job of representing this material. The chemistry of the 
gas plays a significant role in determining the temperature 
only if one looks at colder, denser gas, with much higher 
C and CO fractions, but even here, changing the chemistry 
changes the gas temperatures by at most a factor of two. 
These small temperature changes have little effect on the 
larger-scale gas dynamics, and so the resulting density PDF 
is the same regardless of which chemical model is adopted. 
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APPENDIX A: REVISED THERMAL MODEL 

Our basic treatment of the thermal evolution of the gas is 
the same as in [Glover et al. (20101, and in the interests of 
avoiding duplication, we refer the reader interested in the 
full details of our model to that paper. In this Appendix, 
we restrict ourselves to a discussion of the improvements 
we have made to the basic [Glover et al.| ( |2010| ) model. To 
briefly summarize, we have made three changes. First, we 
have modified our treatment of the adiabatic index 7 and 
the mapping from temperature to internal energy (and vice 
versa) to account for the fact that the internal energy levels 
of H2 are not excited at very low temperatures. Second, we 
have relaxed the assumption made in [Glover et al.| ( |2010[ ) 
and Glover & Mac Low (2011 1 of a constant dust tempera- 



ture, and now solve self-consistently for Tdust, thereby also 
improving the accuracy of our treatment of energy trans- 
fer from the gas to the grains. Third, and finally, we have 
improved our treatment of CO rotational cooling by extend- 
ing the tabulated cooling function used in the code to lower 
temperatures. All three changes are discussed in more detail 
in the sections below. Together, these changes dramatically 
improve the accuracy with which we can model very cold 
gas in our simulations, and allow us to remove the artificial 
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temperature floor at 10 K that was adopted in our previ- 
ous studies. We now should be able to accurately model the 
temperature of the gas down to values as low as 5 K, and 
we find in practice that there is very little gas below this 
temperature in any of our simulations. 

Al Equation of state 

The specific internal energy e for an ideal gas with a temper- 



ature T and a partition function z can be written as ( Boley 
et al |2007| l 



(Al) 



fl ol 

where R = k/rUp is the gas constant, k is Boltzmann's con- 
stant, TUp is the proton mass and fi is the mean molecular 
weight (in units of the proton mass). Because the gas is 
ideal, its specific heat capacity at constant volume, c„, can 
be written as Cv = de/dT. Many previous numerical studies 
of star formation have assumed that c„ is independent of 
temperature, and hence that e = CvT. Unfortunately, this 
assumption is not true for molecular g Boley et al.| 

( 2007 1 have recently pointed out. In very cold molecular gas. 



essentially all of the H2 molecules sit in the ortho or para 
ground state, and the gas behaves as if it were monatomic, 
with Cv = (3/2)_R//i. However, as the temperature of the gas 
is increased, first the rotational and then the vibrational en- 
ergy levels of II2 become populated, and hence c„ changes. 
For this reason, we do not assume a constant in these sim- 
ulations. Instead, we have constructed a set of lookup tables 
that give e as a function of T and XH2 (the fractional abun- 
dance of H2), T as a function of e and , and the adiabatic 
index 7 = Cp/c„ as a function of e (or T) and xu2- When- 
ever we must convert from e to T (or vice versa), or require 
a value for 7, we compute it by interpolating between the 
values stored in these lookup tables. In constructing these 
tables, we have assumed that the H2 ortho-to-para ratio has 
its thermal equilibrium value. 

We have verified that by using this approach, we can 
reproduce the values shown for 7 and e/R in Figures 1 & 2 of 
Boley et al. ( 2007| |. However, we caution that the values for 
these quantities that are used in our actual simulations differ 
slightly from those in these figures, owing to the fact that 
we account for the presence of helium in the gas, whereas 
Boley et al.l (120071) do not. 



A2 Dust cooling and the dust temperature 

We assume, as in most studies, that the dust is in thermal 
equilibrium and solve for the dust temperature Td by finding 
the value that satisfies the equation of thermal balance for 
the dust: 



Ad 



+ r„d + Fh^ = 0. 



(A2) 



Here Text is the dust heating rate per unit volume due to 
absorption of radiation from the interstellar radiation field. 
Adust is the radiative cooling rate of the dust, Fgd is the 
net rate at which energy is transferred from the gas to the 
dust by collisions, and Fna is the dust heating rate per unit 
volume due to H2 formation on the grain surfaces. 

To compute Text , we follow | Goldsmith] ( |2001[ ) and ex- 
press it as the product of a optically thin heating rate, Fext,o, 



and a dimensionless factor, that represents the attenua- 
tion of the interstellar radiation field by dust absorption: 



Foxt — X-TcxtjO- 

The optically thin heating rate is given by 



(A3) 



dv, (A4) 

Jo 

where T> is the dust-to-gas ratio, p is the gas density, is the 
mean specific intensity of the incident interstellar radiation 
field (ISRF), and is the dust opacity in units of cm^ g~^. 
For our ISRF, we used the radiation field given in [Black 
( 1994 1, which is a combination of data from Mathis, Mezger 
& Panagia ( 1983 1 at short wavelengths and Wright et al 



( 1991 1 at long wavelengths. For our dust opacities, we use 
the values given in Ossenkopf & Hcnning ( 1994 ) for grains 
that ha ve not coagulated and that are coated with thick ice 
mantles [Ossenkopf fc Henning ( 1994 1 quote opacities only 



for wavelengths longer than 1/im, and at shorter wavelengths 
we use values taken from , Mathis, Mezger fc Panagia^ (1983] ). 
Finally, for V we take the standard value for solar metallicity 
gas. With these choices, we find that 



Fcxt.o = 5.6 X 10 ^^n erg s ^ cm ^, 



(A5) 



where n is the number density of hydrogen nuclei. 

To compute the attenuation factor x, we solve the equa- 
tion 



x(^h) 



47r J^K^ exp (— K^E) du 
47r J^Ki, du 



(A6) 



where E = 1.4mpA'^H, rnp is the proton mass, and A^'h is the 
number density of hydrogen nuclei. We solve this equation 
for a wide range of different values of A^h , and tabulate the 
results. In our simulations, we can then calculate x using our 
standard six-ray approximation. Recall that this yields, for 
each zone in the simulation, the column density of material 
between the zone and the edge of the simulation volume 
along six rays that we take to be parallel to the coordinate 
axes, for convenience. We can associate a value of x with 
each of these column densities, and the appropriate mean 
value to use in our calculation of the dust temperature is 
then simply the arithmetic mean of these values. By using 
an approach based on lookup tables, we avoid having to 
perform the frequency integration in Equation |A6| during 
the simulation itself, for a great saving in computational 
time. In Figure |A1[ we show how x varies with A'^h for a 
wide range of values of the column density. 
For the dust cooling rate Adust, we solve 



Adust (Td ) =47rOp / B,{Td)K,du. 



(A7) 



where Bv{Td) is the Planck function for a temperature Td, 
and where our values for k,^ were the same as we used to 



* Strictly speaking, we would not expect grains exposed to the 
full interstellar radiation field to be ice-coated, as observationally 
we find that ices appear only for visual extinctions Ay > 3.3 
( [Whittet et al.||1988[ l. In practice, however, our results for Td 
are not particularly sensitive to the particular choice of opacities 
from Ossenkopf &; Hcnning] ( [1994^ , and so we can safely ignore 
this complication. 
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Figure Al. Dust attenuation factor Xi plotted as a function of 
A^H, tlie column density of hydrogen nuclei. 



calculate Text above. We find that the resulting cooling rate 
is well fit by the function 

-3 



Adust (Td) =4.68 X 10 



Td n erg s 



(A8) 



for dust temperatures 5 < Td < 100 K, which more than 
spans the range of values considered in this study. We do 
not account for the absorption by dust of its own thermal 
radiation when computing Adust (Id), as we do not follow the 
evolution of the gas to high enough densities for this effect 
to become important. 

For the net rate of energy transfer from gas to dust, 
Fgd, we adopt the expression ( Hollenbach fc McKee|1989 1 

Tgd = 3.8 X 10 

where 

OL = 1.0 



-33rpl/2 



T^'^aiT - Ti)n^ erg s"^ cm" 



(A9) 



0.8exp(— ) 



(AlO) 



Finally, for the rate at which dust is heated by the for- 
mation of H2 on grain surfaces, Fhs , we have the expression 



Fh2 = 7.2 X 10 /dust-RHa ergs cm 



(All) 



where J?h2 is the formation rate of H2 per unit volume ( Hoi 



lenbach & McKee 19791, and /dust is the fraction of the 



4.48eV H2 binding energy that is absorbed by the grain. 



We assume, following Takahashi & Uehara (20011, that 
/dust = 0.04. 



A3 CO cooling 

We have extended the range of temperatures for which we 

([2OIOI), 



can calculate CO cooling rates. In Glover et al 



made use of the CO cooling function developed by Neufeld fc] 
|Kaufman| ([l993| and [Neufeld, Lepp fc Melnick] |l995| ) and 
hence were limited by the range of their tabulated data, 
which covers the temperature range 10 < T < 2000 K for 
CO rotational cooling. In our present study, we have ex- 
tended their cooling function down to 5 K, using CO de- 



excitation rates taken from Flower ( 2001 1 and Wernli et al. 
( |2006| ) and made available in a convenient form by LAMDA, 



the Leiden atomic and molecular database (|Sch6ier et al 
2005|. 



